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ABSTRACT 

In this work we use the radiation hydrodynamic code TRAMP to perform a two- 
dimensional axially symmetric model of the layered disc. Using this model we follow 
the accumulation of mass in the dead zone due to the radially varying accretion rate. 
We found a new type of instability which causes the dead zone to split into rings. This 
" ring instability" works due to the positive feedback between the thickness of the dead 
zone and the mass accumulation rate. 

We give an analytical description of this instability, taking into account non-zero 
thickness of the dead zone and deviations from the Keplerian rotational velocity. The 
analytical model agrees reasonably well with results of numerical simulations. Finally, 
we speculate about the possible role of the ring instability in protoplanetary discs and 
in the formation of planets. 
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1 INTRODUCTION 

The layered disc model was proposed bv lGammiel il996t) to 
account for accretion-related phenomena in T Tauri stars. 
He assumed that the angular momentum is transported by 
the magne to-rotational instability, commonly referred to as 
the MRI dBalbus fc Hawlevl Il99ll) . However, in the outer 
disc (beyond ~ 0.1 AU) the temperature and the ioniza- 
tion degree is so low that the gas is not well coupled to the 
magnetic field and the MRI decays. There, the only parts 
of the disc in which the MRI can operate are the surface 
layers that are ionized by cosmic rays (ionization due to x- 
ra y quanta emitted by the central star was also considered; 
see iGlasseold. Naiita fc Igeall997l) . Sandwiched between the 
active surface layers is an MRI-free, and, consequently, non- 
viscous area near the mid-plane of the disc, commonly re- 
ferred to as the dead zone. 

An interesting property of layered discs is that, in gen- 
eral, the accretion rate M is a function of the radius r. The 
specific form of this function depends on the mass-weighted 
opacity. However, M increases with r if the opacity does 
not depend on the density and increases with the temper- 
ature not faster than T 2 . This is true in the range 203 K 
< T < 2290p 2 / 49 K, where the opacity is dominated by ice- 
free grains, and also between 167 K and 203 K, where the 
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opacity drops due to the sublimation of ices as T grows. 
At temperatures below 167 K the opacity varies as T 2 , and 
the accretion rate is constant as a function of r (in this pa- 
per we use opac ities of gas and dust mixture according to 
iBell fc Lir]|l994l) . 

When M increases with r, an annulus centred at ro re- 
ceives more mass per unit time from the outer disc (r > ro) 
than that it loses to the inner disc (r < ro). As a result, the 
mass accumulates in the dead zone. Eventually, at one or 
more locations in the dead zone the surface density becomes 
so large that the heat released by the accretion of the accu- 
mulated matter pushes the temperature to a level at which 
the collisional ionization can restore the coupling between 
the gas and the magnetic field. In such case a triggering 
event, e.g. perturbation due to the passage of the companion 
star, heat flux from the in ner active disc or gravitational in - 
stability of the dead zone jArmitaee. Livio fc PringlefeOOlD . 
could start the accretion and "ignite" the dead zone, mak- 
ing at least part of it active. Consequently, the rate at which 
mass is accreted by the central star coul d increase dramati- 
cally. This mechanism was suggested bv|Ganrmie] l^96 j) to 
expla in the FU Ori type outbursts llHartmann fc Kenvonl 
11996ft . 

The layered disc model was further developed bv iHurel 

(2002), who noted that the non-zero thickness of the dead 
zone is an important factor influencing the structure of the 
active surface layers (this is because their structure depends 
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on the vertical component of gravity, which in turn depends 
on the distance from the mid-plane). On the other hand, 
detailed MHD simulations of surf ace layers performed in a 
3D shearing box approximation bv lFleming fc Stond j2003t) 
showed that MRI-driven, non-axisymmetric density waves 
can propagate far into the dead zone. As a result, a purely 
HD turbulence is excited there, providing an effective vis- 
cosity. 

Our original intention was to study how the evolution 
of the dead zone may be affected by the radiative heat- 
ing from the inner active part of the disc. To that end, we 
performed extensive two-dimensional (axisymmetric) simu- 
lations of the evolution of layered discs, allowing for radia- 
tive energy transfer in both radial and vertical directions. 
Unexpectedly, we found that the dead zone is unstable in a 
way that has not been reported before; namely, it tends to 
decompose into rings. In the present paper we illustrate this 
"ring" instability with numerical simulations, and we discuss 
it analytically, providing a physical explanation of the ob- 
served phenomena. The remaining results of our simulations 
will be reported in a forthcoming paper. 

The outline of the present paper is as follows. In §2 
we briefly describe the numerical code and we list the as- 
sumptions underlying our simulations. A layered-disc model 
whose dead zone decomposes into rings is presented in §3. §4 
contains an analytical discussion of the ring instability. Fi- 
nally, in §5 we summarize our results and discuss the effects 
of the ring instability in more sophisticated disc models. 



2 NUMERICAL METHODS 

The simul ations are perform ed with the help of the 3-D code 
TRAMP jKlahr et allll999l) . The equations of continuity 



^ + V.(pv)=0, 
momentum conservation 

^ + (V • V) P V = -VP + pV$ + fccn + V ■ T 

and internal energy 
-8T 



Cvp 



M 



+ (v- V)T 



-PV ■ v + T : (Vv) 



(1) 
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are solved in spherical coordinates (r, 9, (f>) using an explicit 
operator-splitting method. $ is here the gravitational po- 
tential, fccn is the centrifugal force, T is the viscous stress 
tensor, Vv is the rate of strain tensor composed of deriva- 
tives of velocity components, T is the temperature (assumed 
to be the same for gas, dust and radiation), c„ is the specific 
heat, and the colon denotes a double scalar product of two 
tensors. The self-gravity of the disc is neglected, so that 
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where is the mass of the central star. We use the equation 
of state of the ideal gas 



(4) 



The radiation transport is treated at the end of each 
hydrodynamic time-step by solving the equation 

dE r 
dt 



-V • F 



(•») 



where fee is Boltzmann constant, fi is the average molecular 
weight of the disc gas, and ran is the proton mass. 



where E r = aT 4 is the radiation energy density, and F is 
the radiative flux. We adopt 

F = -— VE r , (6) 

where k is the Rosseland mean opacity for the gas-dust 
mixture jBell fc LirJll994h . and A is the flux limiter used 
to int erpolate between optically th in and optically thick 
cases IlLevermore fc Pomraninal98 If) . Equations @ and (|SJ 
are solved implicitly using an iterative SOR method. After 
the convergence has been achieved, the updated tempera- 
ture is calculated from the updated radiation energy den- 
sity. The advection routine for mass, momentum and en- 
ergy is based on a second -order monotoni c transport scheme 
originally introduced by Ivan Leeil jl977f) and optimized by 
iKlev fc Henslerl (I1987T). For more details about the numeri- 
cal code we refer to lKlahr et al.1 dl999h . 

The models are axially symmetric. We allow for the flow 
through the midplane of the disc, i.e. we do not impose a 
reflecting boundary condition there. Grid points are spaced 
logarithmically in r, resulting in a progressively smaller ra- 
dial extent of grid cells near the centre of the disc, where 
their vertical extent is also progressively smaller. Thus, the 
shape of grid cells is more nearly regular throughout the 
grid. 

In the inner active part of the disc, approximate initial 
distributions of density a nd temperature are obtain ed from 
analytical a-disc models dShakura fc Sunavevlll973l) assum- 
ing a constant temperature profile in the direction perpen- 
dicular to the midplane. The outer, lay ered part of the disc is 
initiated with the analytical model of iGammid (|l996). The 
two solutions merge at the radius ruz, where the mid-plane 
temperature of the layered part reaches 1000 K. The accre- 
tion rate of the layered part at ruz defines the accretion 
rate in the inner active part. Similarly, the surface density 
of the dead zone is calculated assuming continuity of the 
total surface density at ruz- Initially, the surface density of 
the dead zone is constant in r. The initial rotational veloc- 
ity is chosen so as to balance the gravity reduced due to the 
radial pressure gradient (in a thin disc the rotation is nearly 
Keplerian) . Other components of velocity are set to zero. 

At the outer edge of the disc mass is injected at every 
time-step at a fixed rate by the following procedure: 

(i) densities and temperatures are copied into the ghost 
zone from the adjacent active cells 

(ii) densities in the ghost zone are normalized to the ini- 
tial surface density 

(iii) a uniform radial velocity is set across the ghost zone 
such that the accretion rate is the same as in the analytical 
model. 

The angular rotational velocity at both inner and outer 
boundary is extrapolated using a r~ 3//2 power-law. At the 
inner edge of the disc, and at the upper and lower boundary 
of the grid, a free outflow boundary condition is imposed. 
A constant temperature of 10 K is maintained at the upper 
and lower boundary. 



At each time-step the location of the dead zone is found 
based on two conditions that have to be fulfilled simultane- 
ously: 

(i) T < Tn m , where Tu m is the minimum temperature at 
which the coupling between the gas and the magnetic field 
still operates 

(ii) E > E a , where E is the column density integrated 
from the surface of the disc along a line perpendicular to 
the midplane, and E a is the column density ionized by the 
cosmic rays. 

Beyond the dead zone the visco sity coefficient is defined ac - 
cording to the a-prescription of lShakura fc Sunavevl lll973l) 

v = ac B H a , (7) 

where a is a dimensionless parameter, c s is the sound speed 
and H a is the thickness of the active layer which is de- 
termined from the vertical distribution of density for each 
radius at each time-step. We chose this particular form of 
Shakura-Sunyaev prescription because it is closer to the an- 
alytical, vertically averaged model. Within the dead zone we 
set v = 0. 

To avoid numerical problems, we introduce a density 
limiter: whenever the density falls below p m in (which may be 
different for different models), it is doubled, and the temper- 
ature is adjusted so as to keep the pressure unchanged. This 
procedure leads to the formation of a low-density "atmo- 
sphere" surrounding the disc. To prevent it from collapsing 
onto the disc, we artificially damp vertical and radial veloci- 
ties in this region. However, this artificial readjustment only 
affects a negligible number of grid cells in the " atmosphere" , 
and the overall evolution of the disc is not influenced. 



3 RESULTS OF SIMULATIONS 

We obtained a broad sample of models with different pa- 
rameters, which will be presented elsewhere. Here we de- 
scribe a representative model with "canonical" parameters 
H = 2.353, Ti im = 1000 K, E a = 100 gem" 2 and a = 0.01. 
Note that at T = 1000 K the ionization degree increases 
by five orders of m agnitude due to ionization of potassium 
iUmebavashil ll983T) . and the magnetic Reynolds number, 
which is directly proportional to the ionization degree, ex- 
ceeds unity. E = 100 g c m" 2 is the standard stopping co lumn 
density of cosmic rays dllmebavashi fc Nakanolll98ll) . and 
0.01 is a value of the vis cosity parameter often adopted fo r 
protoplanetary discs (e.g. lHawlev. Gammie fc Balbusll995l) . 
The simulations were performed on a grid of 256 x 45 points 
in (r,6), with r varying between 0.05 and 0.35 AU and 6 
varying between —5° and +5°. The model was integrated 
for 124 yr, i.e. for 600 orbits at the outer edge of the disc. 

The initial model is far from thermal equilibrium, how- 
ever it quickly relaxes due radiative heat diffusion. The den- 
sity and mass flux across the computational domain shortly 
after the relaxation, at the beginning of the instability and 
at the end of the simulation are shown in Fig. one clearly 
sees how the dead zone successively breaks into more and 
more rings. 

As we have not imposed any explicit perturbations, 
there is a good reason to believe that the rings result from 
a linear instability of layered discs. The mechanism of this 
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v = ac s H a Vj < v 2 




Figure 2. Upper part of the figure shows the layered disc struc- 
ture with the ring-like perturbation. The higher gravitational 
force compresses the elevated part of the active layer, making 
H a i smaller than the unperturbed value H a ^2- Therefore also 
v\ < vi- The radial profiles of average viscosity and accretion 
rate are shown at the bottom of the figure. 



instability is illustrated in Fig. [5] Assume a small axially 
symmetric enhancement of the surface density at a radius 
ro (which, of course, must be accompanied by an increase of 
the dead zone thickness Huz)- Since at higher distances from 
the mid-plane the vertical component of gravity is stronger, 
the active layer (whose column density is at all times fixed 
by cosmic rays) must get thinner at ro. 

In the standard prescription, the viscosity coefficient is 
proportional to the local thickness of the active zone H a (this 
is motivated by the idea that H a is a natural scale that lim- 
its the maximum size of the edies of the MHD turbulence) . 
Therefore, smaller H a results in a lower viscosity in the el- 
evated part of the active layer. The accretion rate, which 
depends on the derivative of viscosity (see Eq. |H| below) , 
decreases at the inner edge of the ring. This causes a bottle- 
neck effect in the accretion flow, and the mass accumulates 
in the ring. On the other hand, the higher accretion rate 
at the outer edge of the ring makes the ring more compact. 
This positive feedback between the dead zone thickness and 
the mass accumulation rate leads to the formation of the 
rings. 

The rings formed in the simulation are shown in the 
radial profile of the surface density in Fig. |3] The accretion 
rate profile exhibits the described minima at the inner edges 
of the rings and maxima at the outer ones. The plot also 
shows the deviations from the Keplerian angular velocity 
- it can be seen that the inner edges of the rings rotate 
at a super-Keplerian velocity. Thus, the rings may capture 
and concentrate the radially drifting boulders of meter-size, 
preventing them f rom accretion onto the central star (e.g. 
iKlahr fc Lmll2000l) . Such conc entration of solids in p e aks o f 
gas density wa s obser ved by lHaghighipour fc Boss] (120031) 
and iRice et all (I2004T) in the case of spiral density waves 
formed by the gravitational instability. 
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Figure 1. A detail of the computational domain around forming rings at the end of the simulation. Contours 
indicate density (the levels are 5 X 1CT 10 , 10~ 9 , 2 X 1CT 9 , 5 X 1CT 9 , 1CT 8 , 2 X 10~ 8 and 5 X 10~ 8 gem" 3 ), 
arrows denote the mass flux (for the sake of lucidity their sizes are proportional to the square root of the 
mass flux). The thick dash-dotted line marks the dead zone. The surface of the disc (i.e. the boundary 
between the disc and the atmosphere) is marked by the dashed line. 
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0.35 



Figure 3. Radial profiles of surface density E (top), the accretion 
rate M (middle) and difference to the Keplerian angular velocity 
SI — Qq (bottom). 



4 ANALYTICAL DESCRIPTION OF THE 
RING INSTABILITY 

4.1 Basic assumptions and definitions 

Let us consider a layered disc consisting of a dead zone with 
thickness 2Huz and two active surface layers with thickness 
H a each. The surface layer accretes at a rate 

M = 6nr 1/2 -?-(2Y, a vr 1/2 ), 



dr 



(8) 



where E a is the column density of the surface layer, v is the 
kinematic viscosity, and r is the cylindrical radial coordi- 
nate. Equation JHJ is a direct conseq uence of angular mo- 
mentum conservation (G el ll996l) . Assuming that the 
accretion energy is radiated locally, we obtain the standard 
formula for the effective temperature T e 



^ a uQ 2 = aT e 4 , (9) 

where a is Stefan-Boltzmann constant. In the optically thick 
approximation the temperature T; at the boundary between 
the dead zone and the active layer is given by the formula 



-rTt 
8 



(10) 



where r = E a /t is the optical depth of the active layer, and 
k is the Rosseland mean opacity llHubenvlll990t) . In gen- 
eral, the opacity coefficient is a complex function of density 
and temperature. In a protoplanetary disc it can be approx- 
imated by a set of power laws 



k(p, T) = Kop a T b 



(11) 



where the constants kq, a and b have diff erent values in dif- 
ferent opacity regimes feell fc Linlll994l) . At temperatures 



lower than 2290p 2//49 K (i.e. nearly everywhere in a proto- 
planetary disc) k practically does not depend on density, so 
we assume a — 0. As before, we assume that the viscosity 
coefficient is given by equation J7J. 

Let 6 be the ratio of the total disc half-thickness to the 
active layer thickness H a : 



8 = 



(12) 



From the equation of hydrostatic equilibrium in the direc- 
tion perpendicular do the midplane of the disc we get an 
approximate formula 

^ = n 2 5H 2 = c 2 si , (13) 
Pi 

where Pi, pi and c s i is, respectively, pressure, density and 
sound speed at the boundary between the dead zone and 
the active layer. Combining 1131 with equations @ and 1101 . 
and using Q with c s = c s i we get 



Ti = 



3ko 9 ks 
8 4<r prnu 



> a 3-bS 2(3-b)f)3-b 



(14) 



Inserting equations @ - (1141 into equation (JSJ, we arrive 
at the following formula for the accretion rate in the active 
layer: 

d 



M = Mr 
where 

M = 12tt 



" <$2(3-f>) r 



1/2 ) , 



3ko 9 



k B 



pniH 



V 3 — b „ T5— 



(15) 



(16) 



In the following, we shall derive an equation which relates 
the angular rotational velocity Q to the thickness of the disc. 

Let us assume that the disc is thin (8H a <C r) , and 
neglect the dependence of SI on z (see e.g. IUrpinlll983l) . We 
may write 



SI 2 = Slo + 



1 OPrr 



Tpm 



Or 



(17) 



where P m = c 2 ; p m , p m and c s i are, respectively, mid-plane 
values of pressure, density, and sound speed. Note that be- 
cause there is no heat generation in the dead zone, the mid- 
plane sound speed is the same as at the boundary between 
the dead zone and the active layer (i.e. at each r the dead 
zone is isothermal along a line perpendicular to the mid- 
plane) . 

The mid-plane density p m is given by the vertical hydro- 
static equilibrium 



p m = Pi exp 



2c 



= pi exp 



-I) 2 
25 



and for i > 1 we have 
'5 
v2 



p m = Pi exp 



(H 



(18) 



(19) 



where it was assumed that the mass accumulated in dead 
zone is already large (i.e. 5 ^> 1), which allows us to neglect 
terms of the second order in 1/8. Inserting p m and P m into 
equation I17t we obtain 



if 



Q 2 , £si_cW 1 dcsj c si dpt 
2r dr r dr pir dr 



(20) 
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The sound speed c s j only weakly depends on disc thickness 
(c 2 . ^ S~ 1/B for 6 = 0.5). The dependence of p t on 8 is 
even weaker (p< ~ <5 1/10 for 6 = 0.5). Therefore, we assume 
that for a small perturbation of 8 the last two rhs terms 
in equation I20H . which are proportional to derivatives of 
Cgj and log(pi), are small compared to the second rhs term, 
which is directly proportional to the derivative of 8. 
The final relation between Q and 8 is 

« 2 =^ + |^. (21) 
2r Or 

Equations 115H and 1211 1 will be used in the next subsections 
to derive the dispersion relation for the perturbations of the 
disc. 



3M -2 + 6-9- 



8nr% 3-6 3-6 



Ah ~ 2 + b - i+b 
O 5=1 A 2 < 3 " b > 



(29) 



and the linearized equation which describes the evolution of 
surface density perturbation with wavenumber k is 

* M r 



2nr 2 



3-2 + 6 -4 + 6 



4 3 - 6 2(3 - 6) 3 -6 



aj-Ah - 2 + b -i"+36 



+ 



In 2 9 _L h 7 — % -8 + 36 -4+b 

■jCsi f + ' n 3-b r 2(3-6) 7^2 

4 3 - 6 2(3 - 6) 
-4 + 6 , 



2(3 - 6) 



-2 + b -10+36 
z — r 2(3-b) ,2 

r Q s J o 



k 1 



(30) 



4.2 The linear analysis 

Let us perturb 8 at a radius rn, assuming that in a small 
region (rn — R, ro + R), R -C r <5 consists of the unperturbed 
part <5n (which in that region may be regarded as indepen- 
dent of r) plus a cosine term with a wavenumber k 



8 — So + S k cos(kR) 



(22) 



where S k is the amplitude of the perturbation. 
According to 1211) . the angular velocity consists of the Ke- 
plerian part Qo and the pressure-correction Sli. Neglecting 
terms of the second order in Oi, and using equation 12 It we 
obtain 

„2 



tt 2 =0,1+ 2Slosli 



fii = 



06 



4rSl dr ' 

The first and second derivatives of 8 and Q axe 

2 

^-8 k k 2 cos(kR) 



(23) 



-8 k k sin(kR) 



:m (l 

2r 



OS 
BR 

- ~8 k k 2 cos(kR) s!n-i5" 



an 

8R 



dR 2 



OR 1 



^ + ^S k k 3 sin(kR) 



(24) 



Since we only want to obtain a rough idea about the growth 
rate of the perturbation with a specific k, we may approx- 
imate cos(kR) with 1 and sin(fci?) with 0. Then the disc 
thickness 8, the angular velocity fl and their derivatives are: 

' rS = -S k k 2 (25) 



08 



8 = 8 + 8 k , -=0 



OR 2 



00 3fi c£ 2 9 2 fi 15S1 

n = n °' ■0R^~2r-^^ 5kk < a^ = -^- (26) 

Due to mass accumulation, the surface density of the dead 
zone Edz grows at a rate 



(27) 



t m = — — 

2nr Or 

Inserting equation 1150 into l|27fl and using approximations 
l|25[l and we obtain 



M 



+ 



+ 



+ 



27rro 

-4 + 6 
2(3 - 6) 
-2 + 6-5 



3 —9 4- 6 - 5 + 2i > -i+b «n 
2 3- 6 9i? 



~ 2 + '' -10 + 3b f) 2 X 



2b 



3-6 2(3-6) 
_2 + 6 



rofl 



-w> -4+6 /an\ 2 



\0Rj 



OR 2 



3-6 

The linearized unperturbed part of the last equation is 



(28) 



The surface density of the dead zone Edz is related to the 
disc thickness 8 through 



= 2 



p(z)dz = 2 



o 



2 / pi exp 
Jo 



2c 2 



Pm exp 



exp 



_D 2 £ 
2e 2 



d« (31) 



where we used the mid-plane density given by 1191 . The 
integral can be evaluated using an error function to yield 



2nc si ( 8 

— ^ exp U 



1 erf 



(32) 



Remembering that 8 3> 1 we get 

E DZ = E a V8 exp (~ - lj 

Differentiating the last equation with respect to time leads 
to a formula connecting Edz with 8 



(33) 



E DZ = E a -f (i + l)exp(£-l) 



(34) 



whose unperturbed part and equation I29H can be combined 
into an equation 

_ _ -7+26 

So 



3M 



2 + 6-9 + 46 ^^^3(3^ 



3-6 



47rr 2 E a 3- 



(35) 



which describes the evolution of the unperturbed disc thick- 
ness due to the accumulation of mass. Next, from equation 
1301 . and equation 1341 1 written for a specific wavenumber 
k we get the dispersion relation which describes the growth 
rate of the perturbation with the wavenumber k. 



8 h = 



M A So 
^El^l^T 



3 _2 + 6 -9 + 46 ^,^ 
43-6 3-6 



<5o + 2 



+ 



-4 + 6 1 



+ 



2(<5 + l) 2(3- 6) 8 

34 -2 + 6 7 - 36 - 
4 3- 6 2(3 - 6) 

_4 + 6 - =2±* 



8 + 36 -7+26 

^S^k 2 



2(3 - 6) 



(36) 
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Figure 4. Growth rates of the four inner-most rings depending 
on the unperturbed disc thickness So in comparison to growth 
rates calculated from equation 1361 . Growth rates in the numeri- 
cal simulation were determined by measuring 8 and <5 (u> = S/S). 



The perturbation growth rate diverges at short wavelengths. 
This is because our simple analysis does not include any 
damping mechanisms. In reality, however, thin rings would 
diffuse due to the thermal motion of particles. 

Note also that that the radial pressure scale height can- 
not become smaller than the vertical one. Therefore, rings 
with circular r-z profiles should form, which is in agreement 
with our numerical model (see Fig. 0. 

Fig. [I] shows growth rates of the disc thickness 8 in 
the four inner-most rings. They are compared to the growth 
rates given by the dispersion relation 1361 . The growth rates 
in the analytical model are more than one order of magni- 
tude higher. The possible explanation of this discrepancy 
can be that the analytical model does not allow for heat 
transfer from the ring into the active layer above it or for 
convective flows that mix the mass inside the rings. The 
radiative transfer in the radial direction may also be im- 
portant. All those processes make the dip in the average 
viscosity shallower and effectively decrease the growth rate 
of the instability. 



4.3 The influence of stellar irradiation 

In this section we estimate how strongly the ring instability 
is affected by the radiation flux from the star. In the preced- 
ing section the stellar flux was not included self-consistently 
since the linear analysis presented there serves as the com- 
parison model for the numerical simulation which does not 
include irradiation. On the other hand, with the irradiation 
included explicitly the calculations would become too com- 
plex (or just impossible) to perform. 

Firstly, we estimate how important the irradiation is 
for the unperturbed disc. We do it by computing the change 
of the temperature at the boundary between the dead zone 
and the surface layer (from which the thickness of the surface 
layer and the viscosity are determined). 

Stellar heating can be included in Eq. 1101 in the fol- 
lowing way 



where T e is the effective temperature resulting solely from 
the viscous dissipation, is the temperature of the stellar 
photosphere, and W is the dilution factor which depends on 
stell ar radius, dista nce from the star, and geometry of the 
disc jHubenvlll990ll . For a conical disc (in which the aspect 
ratio H/r does not depend on r) and r > i? t we have 



W = 



37t ( r ) 



(38) 



(see e.g. lRuden fc Pollacklll99lh . 

Evaluating the rhs terms in Eq. 1371 for our disc model 
and typical parameters of a T Tauri star (7?* = 3 Rq, 
= 4400 K) we find that the first (viscous) term always 
dominates. More susceptible to effects of irradiation are flar- 
ing discs, in which H/r ~ r 1 . The flaring index 7 depends on 
the model, e.g. a vertically isothermal model in which the 
intercepted stellar flux is equal to the blackbody emissio n 
from the disc yields 7 = 2/7 JChiang fc Goldreichl Il997h . 
However, even in such model the irradiation term dominates 
only for (r > 10 AU), where the grazing angle becomes suf- 
ficiently large. 

Effects of irradiation may be more important for the 
ring instability itself. The inner edge of a growing ring is 
more strongly heated by stellar radiation, because the graz- 
ing angle a gr is larger there. As a result, the temperature 
and the thickness of the active layer increase locally, leading 
to an increased viscosity. Then, the dip in the accretion rate, 
which is responsible for the ring growth, becomes shallower 
or it may even disappear entirely. 

The importance of this effect may be roughly estimated 
by comparing the amplitude of viscosity increase due to stel- 
lar irradiation to the amplitude of viscosity decrease due to 
the stronger vertical component of the gravitational force 
(the latter effect is explained in Fig. 

Let us assume a ring-like perturbation of the disc thick- 
ness with an amplitude Sk and a wavelength A = 2n/k. 
Since the most unstable wavelength is A max ~ So + Sk, let 
us parametrize the perturbation wavelength by the dimen- 
sionless value I = A/A max . The grazing angle at which the 
stellar radiation strikes the inner edge of the ring is 



Sk 2 R* 

l(8o + 8k) 37r r 



(39) 



Evaluating the dilution factor W = Q gr (J?*/r) 2 , inserting it 
into Eq. 13 71 . and combining it with Eqs. 1131 and 
we obtain the viscosity in a form 



v(8k,8o,r,l) = iVi(r)(<5o + <5fc) 2(3 - b) + N 2 (r) 

8k 



x (80 + 6k) 



-1/2 



l(S + Sk) 



+ N 3 (r) 



1/4 



(40) 



where 
Ni 



3Ko 9 \ 3^6 £-b / fc £ 

1 a 3 - b 



if = gTT e 4 + WTt 



(37) 



N 2 (r 
and 
N 3 {r 



(Sua 9 

a k B /R^- 1 ' 2 



\ /im H 



4—b 



-2 + b 

fi^ir t (41) 



f2 (IMh 

37t r 



(42) 



(43) 
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Figure 5. The contours show the critical amplitude 5fe jC rit(5o,r) 
for which v(<5j. cr i t ) = v(0). Beyond their endpoints the solution 
becomes unphysical (5^ > <5oi i.e. the perturbed disk has regions 
with negative thickness). 



Initially, the function u(8k) always grows as 8k increases. 
However, depending on the remaining parameters (So, r and 
I), it may start to decrease and quickly drop below the initial 
value v(8k — 0). Therefore, a very small perturbation of the 
disc thickness is always stabilized by the stellar irradiation, 
but if the amplitude of the perturbation reaches some value, 
the ring instability may start to work. This critical value 
5fc,crit (v(8k,ciit) = v(0)) strongly depends on I: it is smaller 
for higher I, where the grazing angle is smaller. The depen- 
dence on the remaining parameters (r and So) for I = 3 is 
shown by Fig. 

We see that the stabilizing effect of the irradiation may 
become unimportant for broad rings (large I), small So (less 
mass accumulated in the dead zone) and/or at small radii. 



5 DISCUSSION 

We described a new instability in layered accretion discs. 
The accretion rate in such discs is in general a function of 
radius. As a result, mass may accumulate in the non-viscous 
dead zone near the mid-plane of the disc. However, a small 
ring-like perturbation of the dead zone thickness leads to 
the deviation of the rotational velocity which results in non- 
uniform accumulation rate of the mass supporting growing 
of that ring perturbation. This ring instability may eventu- 
ally lead to a decomposition of the dead zone into rings. 

We observed the formation of such rings in the 2D ax- 
ially symmetric radiation-hydrodynamic simulation of the 
layered disc. To illustrate how the instability works, we 
give its approximate analytical description. According to the 
analytical results, the narrowest rings grow most rapidly. 
Therefore, we may expect the formation of the radially thin 
rings whose radial extent will be given just by the thermal 
pressure of the gas. The comparison of ring sizes shows a 
reasonable agreement between the numerical simulation and 
the analytical model. 

The irradiation by the central star may substantially 
decelerate or even stop the growth of the instability in some 
regions of the disc. However, broad rings are less affected, 
and the effects of irradiation become less important in thin 
discs and/or at small distances from the star. On the other 



hand, once the innermost ring has developed, and the flaring 
index 7 is not too high, the disc at larger radii is shadowed 
and more rings may develop there. If the disc is truncated 
(e.g. by the magnetosphere of the star) the shadowing effect 
may also be caused by its inner edge, allowing for an efficient 
growth of the instability. 

The rings created by this instability may be important 
in terms of planet formation, because they can be the places 
where the solid particles (gravel and boulders) accumulate. 
It is because of the higher rotational velocity at the inner 
edge of the ring and the lower one at the outer edge. The 
drag force which makes the grains to move inwards is smaller 
at the inner edge of the ring and higher at the outer edge. 
The solid particles may merge into larger bodies (planetesi- 
mals) necessary for the formation of planets. 

Massive rings are subject to a hydrodynamical insta- 
bility in 3D, e.g. with respect to non-axis ymmetric pertur- 



bations JPapaloizou fc Pringlelll984ll985D . This instability 
occurs if the rotational angular velocity decreases with ra- 
dius steeper than r _v ^. In our simulation, this occurs for 
the innermost ring at a time of ~ 280 yr. In the non- 
viscous case such ring would most likely decay into so-called 
"planets", i.e. large scale vo rtices, as found in numerical 
simulations bv | H awley| Jl987l), a nd fu rther investigated by 
iGoodman. Narvan fc Goldreichl il987f) . The fate of such a 
ring in 3D viscous hydro-simulations including the effects of 
layered accretion still have to be performed. 

Our model assumes zero viscosity i n the dead zone. 
Howe ver, there are some indications iFleming fc Stonel 
2003) that even in the dead one there might be some tur- 
bulence present, induced by waves originating in the active 
parts of the disc. In such a case, the excess of surface den- 
sity could result in a higher accretion rate in the rings, and 
the growth of the rings would be stopped or even they might 
not form at all. This issue will be addressed in a forthcoming 
paper. 
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